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Abstract 



The ability to quantitatively assess the health of an ecosystem is often of great interest to 
those tasked with monitoring and conserving ecosystems. For decades, research in this 
area has relied upon multimetric indices of various forms. Although indices may be num- 
bers, many are constructed based on procedures that are highly qualitative in nature, thus 
limiting the quantitative rigour of the practical interpretations made from these indices. 
The statistical modelling approach to construct the latent health factor index (LHFI) was 
recently developed to express ecological data, collected to construct conventional multi- 
metric health indices, in a rigorous quantitative model that integrates qualitative features 
of ecosystem health and preconceived ecological relationships among such features. This 
hierarchical modelling approach allows (a) statistical inference of health for observed sites 
and (b) prediction of health for unobserved sites, all accompanied by formal uncertainty 
statements. Thus far, the LHFI approach has been demonstrated and validated on fresh- 
water ecosystems. The goal of this paper is to adapt this approach to modelling estuarine 
ecosystem health, particularly that of the previously unassessed system in Richibucto in 
New Brunswick, Canada. Field data correspond to biotic health metrics that constitute the 
AZTI marine biotic index (AMBI) and abiotic predictors preconceived to influence biota. 
We also briefly discuss related LHFI research involving additional metrics that form the in- 
faunal trophic index (ITI). Our paper is the first to construct a scientifically sensible model 
to rigorously identify the collective explanatory capacity of salinity, distance downstream, 
channel depth, and silt-clay content — all regarded a priori as qualitatively important abi- 
otic drivers — towards site health in the Richibucto ecosystem. 

Keywords: AMBI; Bayesian statistics; hierarchical modelling; infaunal trophic index; Markov 
chain Monte Carlo; statistical inference 



1 Existing Methods to Quantify Ecosystem Health 



Assessment of the "health'' of an ecosystem is often of great importance to those 
interested in the monitoring and conservation of ecosystems. Health is a com- 
plex concept often involving many diverse factors, and therefore is not straight- 
forward to quantify. A popular method to estimate the health of an ecosystem is 
through one or more multimetric indices, each of which is a scalar number col- 
lapsed from several indicator variables of health, or metrics. Ecosystem health 
metrics are frequently measures of faunal abundance and diversity. For aquatic 
ecosystems, these biotic metrics often focus on benthic populations — organisms 
living on or in the sediment at the bottom of a body of water — since they are 
useful indicators of underlying health conditions (B ilyard , 1987; Dauer[|1993l ). For 
example, the AZT^ marine biotic index (AMBI) (Borja et a/.[|2000| ) is a quantitative 
measure of health for an estuarine ecosystem based on the sample counts of cate- 
gorised benthos. Its popularity is evident from its use across the globe, including 
Africa ([ Bazairi et al\ |2005[), Asia (|Cai et al\ |2012[ ), Europ e ([Medeiros efal] |2012[ ), 
North America ( [Teixeira et a/.[ |2012), and South America (Mun iz et a/.[|2012 ). 

AMBI and other common multimetric indices, e.g. infaunal trophic index (ITI) 



^Marine and Food Technological Centre ( [http : / /www . azti . es| ). 
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( Word[|1980[), estua rine biotic integrity index ( [Deegan et al.\ 1997), b enthic response 
index ( Smith et~al\ |2001|), b enthic quality index ( [Rosenberg et al\ [2004[ ), infaunal 
quality index ( |Mackie[ 2009[ )^ have the main appeal that they are conceptually sim- 
ple and thus easily interpretable. They also contain a high amount of biological 
content from subject-matter scientists being involved at all stages of the design of 
the index. On the other hand, many of these stages during index construction can 
involve a non-trivial amount of arbitrariness. Consequently, rigorous evaluation 
of index reliability and other quantitative aspects is difficult with conventional in- 
dices: for example, detecting relationships between health and environmental or 
impact-related covariates such as water depth or urbanisation; and formally as- 
sessing the uncertainty in these estimates of health. Recent multi-step approaches 
towards addressing such concerns (e.g. Smith et ah , 2001 ; Johnston et al.\ 2009| ) do 
not address propagation of uncertainty from one step to another, thereby resulting 
in inference that is less reliable than that from an integrated statistical method- 
ology. Chiu & Guttorp ( 2006| ) proposed the SHIPSL approach, a statistically en- 
hanced multimetric index construction scheme that improves various quantitative 
aspects of conventional indices ( Dobbie & Dail[ to appear| , although it and oth- 
ers share unresolved issues such as non-transferability in space or time, and the 
need for follow-up analyses to determine its relationship with non-faunal (abiotic) 
variables in method evaluation or policy-making contexts. 

Recently Chiu et ah ( 2011[ ) devised the latent health factor index (LHFI), a novel 
statistical model-based ecological index aimed to retain the advantages of con- 
ventional multimetric indices while addressing some of their shortcomings. The 
LHFI methodology involves a multi-level analysis of covariance generalised linear 
mixed-effects (regression) model (e.g. [Hoff. 2009) , or ANOCOVA GLMM: instead 
of being treated as measures of health, metrics are regarded as indicators of under- 
lying health conditions; thus, these indicators are regressed as response variables 
upon a latent health quantity (latent since it is not directly observable) which is 
site-specific, forming the main level of the regression; health in turn can be re- 
gressed upon available covariates, such as environmental (e.g. salinity, silt-clay 
content) and impact related (e.g. urbanisation) variables, forming the optional sub- 
level in the model hierarchy. 

With data on metrics and covariates, latent health can be estimated as a scalar, 
so that interpretability is retained; the estimated quantity is the value of the index. 
Additionally, the effect of the covariates on health can be estimated in a single in- 
tegrated statistical framework. Importantly, statistical modelling is what directly 
produces the health index under the integrated LHFI framework, as opposed to 
being employed merely to select relevant metrics before index construction (e.g. in 
Deegan et ah\ [T997| ) or to evaluate the resulting index (e.g. in Borja et ah. 2009[ ). 



Thus, the LHFI is much more rigorous than conventional indices, as its definition 
utilises universal modelling practices for the definition of the index; its hierarchical 
modelling framework also allows for comprehensive statistical inference without 
the need for sequential analyses through which the propagation of uncertainty is 
lost from one analysis to the next. As well, the LHFI model, once specified for a 
set of existing sites, allows for cost effective yet rigorous interpolation of health for 
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a new site: prediction can be accomplished simply with information on covariates 
at this new site, thus bypassing the expensive benthic taxonomic laboratory pro- 
cedures that are required to gather the metric data as required by conventional in- 
dices. These desirable properties are gained without sacrificing biological integrity 
which can be embedded through subject-matter expertise in the identification of 
useful metrics and covariates for constructing the LHFI. It is also straightforward 
to use the LHFI framework to handle data that have certain spatial and/ or tempo- 
ral features, thus resolving the non-transfer ability issue of previous indices. 

Recently, [Schliep and Hoeting ( |2012[ ) integrated formal point-referenced spa- 
tial modelling Banerjee et al. (|2004]) with LHFI principles to model the hierarchi- 
cal relationship among four levels of quantities: five-point ordinal health metrics 
(scaled from "poor'' to "excellent''), latent continuous quantities that determine 
the ordinal metrics, latent health, and drivers of health that pertain to geographi- 
cal and environmental characteristics. They illustrate the type of unified statistical 
inference that can be drawn from such an LHFI-based approach for assessing bi- 
otic integrity of river basins in Colorado, USA. In contrast, Chiu et al, ( 2011| ) and 
Wu| ( unpublished! ) directly model the quantitative health indicators based upon 



which ordinal metrics such as those of [Schliep and Hoeting ( |2012[ ) are defined. 
This avoids loss of information due to the mapping of quantitative health metrics 
to a coarse ordinal scale. 

More information on popular indices, the LHFI, and their properties are given 
by |Chiu gfaT] ( |2011[ ) and | Wu| ( [unpublished| ) . 



2 Previous LHFIs for the Richibucto Estuary 



The LHFI modelling methodology was first demonstrated and validated on fresh- 
water ecosystems by Chiu et al. ( 2011[ ). Following this, | Wu (unpublished) applied 
the methodology successfully to an estuarine ecosystem, utilising the dataset from 
Lu et al\ ( |2008[ ) on the previously unassessed Richibucto estuary in the Canadian 
province of New Brunswick. The data were collected in the estuary at 18 sites 
(Fig. [1|). Sites 2, 4-7, and 14 were closest to active oyster farms. Oyster farming 
activity is perceived to impact site health through its direct influence on sediment 
properties, although Lu et al, ( 2008[ ) report that different biotic indicators show dif- 
ferent types of association with proximity to oyster farms. For example, macroben- 
thic faunal abundance is medium for 5 of these 6 sites but high for Site 2; Shannon's 
diversity is relatively even among all 18 sites with a slight upward trend as the site 
moves away from the upper channel instead of from an oyster farm. This lack of 
obvious association remains despite these authors' efforts to consider certain indi- 
vidual species as separate health indicators. This motivated us to develop LHFI 
models for Richibucto based on indicator metrics (Table [l]) used to construct the 
AMBI and ITI, two popular estuarine ecosystem health indices. Specifically, AMBI 
and ITI metrics are better tailored to estuarine ecosystems than the generic indi- 
cators of abundance, richness, and diversity; and they are more comprehensive 
than indicators based on single species. However, biotic health indicators alone do 
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Figure 1: Map of Richibucto estuary with the 18 monitored sites (labelled 'G' fol- 
lowed by the site number); the straight lines illustrate the calculation of distance 
downstream (DD) for Sites 3 and 5. 



not explicitly reveal the collective impact on overall health from abiotic variables, 
including sediment properties (which, in this case, can be affected by oyster farm- 
ing activity), oceanographic properties (salinity, water temperature, etc.), and their 
interactions. The work by |Wu ( |unpublished| ) is the precursor to our current paper. 
Wu ( unpublished! ) developed two sets of LHFI models, the first using only met- 



rics from AMBI (denoted by LHFI-A), and the second using metrics from both 
AMBI and ITI (denoted by LHFI-A-I). These models were considered in a Bayesian 
statistical framework and implemented using Markov chain Monte Carlo (MCMC) 
techniques. Both sets of models were successful in that they were able to make rea- 
sonable distinctions between health levels at different sites, while allowing rigor- 
ous assessment of reliability. As well, the resulting health estimates were justifiable 
by subject-matter expertise. 



2.1 What May Influence Richibucto's Health? 

AMBI metrics, which for the most part pertain to organic enrichment, were the 
main focus when Wu ( unpublished! ) constructed the Richibucto LHFIs. This was 



^The fraction of organisms with the specified characteristics out of all benthic organisms in the 
grab sample. 



^Neither clearly positive nor clearly n egative. 
^As described by Cromey et ah (2002K 



Table 1: Metrics, based on definitions of the AMBI and ITI, that are used here and 
by |Wu| ( unpublished| ) to construct LHFIs for the Richibucto estuary. 



Preconceived 
association 
w/ health 



AMBI abundanceP metric 



species (including specialist carnivores and some deposit-feeding 
tubicolous polychaetes) very sensitive to organic enrichment and 
present under unpolluted conditions 



species (including suspension feeders, less selective carnivores and 
scavengers) indifferent to enrichment, always present in low densities 
with non-significant variations with time 

species tolerant to excess organic matter enrichment (including 
surface deposit-feeding species, e.g. tubicolous spionids) 

second-order oppotunistic species; mainly small-sized polychaetes: 
subsurface deposit-feeders, e.g. cirratulids 

first-order opportunistic species: deposit-feeders, which proliferate 
in reduced sediments 



+ 



ITI abundance metricP 



1 suspension feeders: feed on detritus from the water column and + 
usually lack sediment grains in their stomach contents 

2 interface / surface detrital feeders: obtain the same types of food as + 
suspension feeders but usually from the upper 0.5 cm of the sediment 

3 deposit feeders: invertebrates (including carnivores); generally feed ± 
from the top few cm of the sediment and feed on encrusted mineral 
aggregates, deposit particles or biological remains 

4 specialised environment feeders: mobile burrowers that feed on — 
deposited organic material; all adapted to live in highly anaerobic 

sediment 



because prior to Wu's analyses, it had been perceived that benthic fauna in Richibucto 
were related to organic enrichment, as well as freshwater input (salinity gradient), 
variability of particle size and topography (channel and water depth) ( |Lu et al. 



2008| ). These non-enrichment aspects of the estuary were observed as covariates 
alongside benthic fauna; see Lu et al. (|2008j) for details. 

As such, in addition to assessing the health of the Richibucto system, the above 
LHFI models were used to investigate which and how covariates may influence 
health as reflected by biotic metrics. As discussed by |Chiu et al. ( 2011[ ) and Wu 



( |unpublished| ), a thorough understanding of the relationship between covariates 
and health is key to rigorous yet cost effective interpolation of site health, and 
could prove to be an enormous asset to scientists. To this end, each of the above 
LHFI-A and LHFI-A-I models had been implemented with a different combination 
of covariates. Several covariates and interactions were found to have a statistically 
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Table 2: Sample correlation coefficients among covariates for latent health of 
Richibucto sites. 





DD salinity 


log(depth) 


log(SC) 


DD 

salinity 
log(depth) 
log(SC) 


1 0.88 
1 


0.16 
0.23 
1 


-0.47 
-0.33 
-0.41 
1 



significant relationship with health]^ For LHFI-A, a model with site-associated co- 
variates log(SC), log(depth), salinity, and interaction log(SC)xlog(depth), and another 
model with a single covariate distance downstream (DD), were the two best-fitting 
models among those investigated. For LHFI-A-I, there were three best models: 
two corresponded to the same sets of covariates as those for LHFI-A, and another 
model with covariates log(depth), log(SC) and their interaction. SC denotes the 
fraction of silt-clay (grains of size <63 /xm) out of the sediment pooled from all (2 
to 3) grab sample replicates. Depth is the distance (m) from the water surface to 
the estuary bed at the location of the site from which grab samples were obtained. 
Salinity (parts per thousand) was measured based on one in situ water sample ob- 
tained at the site. To determine DD (km), first a straight line was extended from 
the western-most site (Site 1) to the eastern-most site (Site 18); DD of any site is de- 
fined as the distance between the site's perpendicular projection onto the straight 
line and Site 1 (Fig.jl]). Other site-specific covariates also considered but found to 
be insignificant or confounded with others were water temperature (^C), time of year 
(September or October), median grain size of sediment, sorting (a unitless measure 
of variability of grain size) and organic content (%). 

However, attempts by |Wu ( [unpublished l to include covariates from various 



best-fitting models together in a single LHFI-A or LHFI-A-I model were unsatis- 
factory; in such combined models, DD remained highly significant, while all other 
covariates and their interactions were no longer significant at a reasonable credible 
level. Indeed, salinity and DD are highly correlated (Table[2]), and it was unsurpris- 
ing that the two are not simultaneously significant. However, no strong correlation 
exist among log(depth), log(SC), and DD (Table |2]), and so why did DD "trump'' 
all others in a combined model, despite non-DD covariates being significant when 
DD was absent? This question had yet to be addressed. As well, while relation- 
ships between health and covariates were quite strong for the LHFI-A models, they 
were less clear for LHFI-A-I models (significance at credible levels ?^60-85% in the 
best-fitting LHFI-A-I models as opposed to >90% for LHFI-A). This indicated that 
the extra data from ITI metrics weakened the overall relationship between health 
and covariates. 

One possible explanation for this phenomenon is that the LHFI construct was 
appropriate for describing health using AMBI metrics and the available covari- 
ates, but ITI metrics have weak ecological relevance to Richibucto. This is plausi- 



^ Significance here refers to the regression coefficient having an interval with a reasonably high 
Bayesian credible level (e.g. 0.8) while excluding the value 0. 
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ble from a qualitative perspective, in light of our prior beliefs as stated in the start 
of Section |2.1| To quantitatively address this, one may determine additional co- 
variates that can be more appropriately paired with ITI metrics, then model these 
alongside the original covariates. However, this would require further field activi- 
ties, and we do not pursue it in this article due to cost constraints. 

On the other hand, it is also possible that a) the ITI data were too noisy for Wu's 
specific LHFI models to detect any patterns in their relationship with the covari- 
ates, or that b) the data were not too noisy, but these LHFI models were inadequate 
for revealing a clear relationship among latent health and covariates. Scenario a) is 
certainly conceivable given the type of study at hand, in which data often involve 
substantial measurement error. If this were the case, there is unlikely much room 
for the proposed models to be improved upon within the same modelling frame- 
work, given that they are already quite ecologically informative. This leaves us 
with b) to consider; indeed, Wu's models were merely preliminary models under 
the general LHFI construct. This could also explain the "trumping'' phenomenon 
in the simpler LHFI models. Specifically, distance likely contained much less mea- 
surement error than the other covariates, it being easier to measure with precision 
than the environmental covariates which are intrinsincally more variable in nature. 
With an LHFI model that was perhaps too simplistic, an effect on health from dis- 
tance could therefore manifest itself more clearly than effects from other covariates 
even if all of them are equally important qualitatively. 

Existing data can be used immediately to address b), but they require an im- 
proved quantitative framework. In light of our prior beliefs, and the fact that the 
environmental covariates contained specialised information that distance did not, 
a more sophisticated LHFI modelling framework yet may prove to be helpful in 
providing a common thread through health, DD and the other ecologically rele- 
vant covariates. 

Thus, in the rest of this article, we focus on addressing the "trumping'' phe- 
nonmenon in the context of b). To do so, we wish to determine if implementing 
either or both of the following will allow us to properly identify the nature of the 
relationship among latent health and the available covariates: 

1. Introduce a covariance structure for the metric effects, instead of indepen- 
dence which was assumed for the preliminary models to reduce computa- 
tional burden. 

2. Introduce additional level(s) to the regression hierarchy based upon the known 
associations between the available covariates. 

Since these steps pertain to different parts of the LHFI model, below we treat each 
as a stand-alone investigation. Moreover, we consider LHFI-A models only: intro- 
ducing extra model complexity for LHFI-A-I models can be impractical for proper 
inference (e.g. via MCMC), given that AMBI and ITI metrics are dependent accord- 
ing to their definitions, so that extra model parameters are required to account for 
this. (For example, see Wu| ( |unpublished| ) for the extra parameters involved even 



when such dependence was only informally accounted for.) 
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3 Modelling Covariance of Metric Effects for LHFI-A 



Without additional data, estimating an unstructured covariance matrix may be im- 
practical due to the considerable increase in new parameters in the model. Instead, 
assuming a structured covariance matrix would be preferred. The form of the LHFI 
model, for instance, might inherently imply a certain covariance structure. Thus, 
we examine what covariance structure might be appropriate for the metric effects 
under previous LHFI-A models, and would then implement them to determine if 
the extra model complexity can bring about more detectable relationships between 
latent health and covariates. 



3.1 Extending the LHFI-A Models for Richibucto 



We follow the notation of Wu ( unpublished[ ). AMBI metrics, denoted by Y in the 



LHFI framework, are abundances of five disjoint taxonomic groups. Due to the 
different preconceived directions of their association with health (Table [l]), we split 
the metrics into two groups: s='—'' for Metrics 3-5 (negatively related to health), 
and s='V for the remaining metrics. In the LHFI model, each member of Group 
s is modelled as a multinomial random variable. The link function for the GLMM 
is a generalised logit for s='-\-'', and an inverted generalised logit for s=''—'\ so 
that large metric values for s='V and reflect good/neutral and poor health, 
respectively. 

More precisely, let lixj(s)xfc/ written Yijks for simplicity, denote the value of the 
jth metric (nested within the sth metric group) for the kth replicate grab sample at 
the ith site. Let N^k denote the cardinality (total number of benthic organisms) of 
the fcth replicate sample at the ith site; and pijs denote the unknown probability of 
an organism from the ith site belonging to the j{s)\\\ taxonomic group. Thus, we 
have multinomial distributions 

- Multinomial(A^^fc;pa+,Pi2+, 1 - Pii+ - Pi2+) , 

\Xi?>k- 5 YiAk- ) Yi5k- 5 Nik — YiSk- — ^lAk- — ^i^k- } \ ^i/c , Pi3- , PiA- , Pi5- 

^ Multmomml{Nik] Pi3-,Pi4-,Pi5-^ 1 - Pis- - PiA- - Pi5-) • 

Next, let Hi denote the latent health of the ith site; and 9s and f3j(^s) respectively 
denote the metric group effect and individual metric effect (both unknown) in the 
regression model. Then, the linear predictor in the LHFI framework is 

ly^J+ = log- ^ , i = l,2, (3) 

1 - Pii+ - Pi2+ 

Pij- = log , j =3,4,5, (4) 

Pij- 

Uijs = Hi + 9s + f3j(s) , s = +,-. (5) 

For (|5]), we model 9s as a fixed effect and take 6'+=0 (as is customary when con- 
sidering one of the categories as baseline) to ensure model identifiability, and we 
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model site health and metric effects as random. Specifically, the latent regression 
of Hi is 



Hi = ao + cx'xi + Si 



Si a 



2 iid 
H ^ 



Normal(0, aj^) 



(6) 
(7) 



where Xi is the vector of a given combination of the aforementioned covariates]^ 
ao and cx are the unknown coefficients of the corresponding latent regression, and 
£i is the normally distributed regression error with unknown variance aj^. 

Note that there is overlap and thus dependency between the two multinomi- 
als of ^ and |Wu ( |unpublished[ ) explained that this dependency is crudely 
accounted for by Os) similarly, the mean-zero /3j(^s) crudely accounts for the depen- 
dency among us within group s. Thus, independent /5j(5)S were assumed. More 
rigorously, we now replace independence of metric effects by 



-MVN(0,S), f3 



(8) 



where y9+ = /52(+)]^ f^- = [/53(_), ;54(_), /35(_)]^ S is the unknown covariance 

matrix for /3, and "MVN'' denotes the multivariate normal distribution. Thus, S+ 
(2x2), S_ (3x3), and S± (2x3) denote the covariance matrices for metric groups 
positively and negatively related to health, and their cross-covariance matrix, re- 
spectively. 

We use relatively diffus^ distributions (with the same parametrisation as Wu 



( unpublished! )) as priors for ao,a.,9- (univariate Gaussian with mean and vari 



ance 100) and aj^ (inverse-Gamma with unit shape and scale). To complete the 
Bayesian modelling hierarchy, we must specify the structure for S, as we now ex- 
plore. 



3.2 Dependence Structures of /5s and us 

The most general form of S is to take it as fully unstructured, and thus it can be 
assumed to have an inverse- Wishart (IW) distribution 

S - IW5 (9) 

where ''IW/' denotes the inverse of a x Wishart matrix with d degrees of 
freedom and scale matrix equal to the identity, which is a relatively diffuse prior 

""in practice, covariate transformation might be necessary to satisfy the linearity of Covari- 
ates (possibly transformed) are then centred to reduce dependence among the as. For a given co- 
variate that is not an interaction, centred data are produced by subtracting from the raw covariate 
data a constant that is (approximately) equal to the observed covariate mean (averaged over i). The 
"centred interaction'' between two covariates is taken to be the product of two centred covariates. 

^Diffuseness of priors reflects the fact that in the absence of data, we have no clear perception of 
the properties of the corresponding unknown quantities. In general, diffuseness reduces the need 
for justification of prior distributional assumptions. 
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for a. d X d unstructured covariance matrix. Then, one can take advantage of ex- 
isting MCMC software such as OpenBUGS ( http : //www, openbugs . inf op for 
straightforward implementation of the LHFI model, although in our experience. 



non-trivial hierarchical centring is essential to improve MCMC mixing (see Chiu 
et a/.[|201l||Wu[ [unpublished, for details) 



However, one concern for assuming ^ is that with a small dataset from 18 
sites each with only 2 to 3 replicate grab samples, an unstructured S m ay be only 
weakly identifiable]^ This issue was encountered by Chiu et ah ( |201l l) for a fresh- 



water benthic dataset that also involved 18 sites with 3 replicates per site, although 
there were nine metrics altogether. To address this concern, we instead consider a 
structure for S which involves fewer unknown parameters. 

To this end, let us momentarily consider a frequentist's viewpoint. 

Without loss of generality, we may drop the subscript s from p and u in (|3])-(|5]) 
since the value of j determines s. Then, given site z, let S"^ denote the 5x5 covari- 
ance matrix whose (j, /)th element is SJ^/=Cov(z/^j, Uijf) which does not depend on 
i. In the frequentist context, one can show that /3 has covariance matrix 

S = - 4 J55 (10) 

where Jdd' isadx d^ matrix of Is. Furthermore, we partition S"" into and 
"±'' blocks accordingly. 

Thus, as an alternative to (|9]), one may wish to consider the structure of S"^ 
while limiting the complexity of the resulting structure ultimately assumed for 
S. For this, we assume that the dependence between the vectors of u^s and u_s is 
adequately addressed by Og in (|5]), and consequently is a matrix of Os. However, 
we allow and S"! to be unstructured. 

The form of ( p^ in light of the above consideration suggests that a reasonable 
structure for S may be 

S+l^ - IW2 + <^J22 , - IW3 + c^Jss , ^±\^ = , (11) 

^ N(0, 100) subject to S being positive definite. (12) 

Our choice of distributions and hyperparameters in ([TT]) and (12) yields relatively 
diffuse priors; it also allows S to be decomposed as the sum of the constant matrix 
J55 and a random block-diagonal matrix whose blocks are unstructured, thus 
giving S a general form that mimics that in ([TO]) while keeping S positive definite. 



Altogether, our extended model comprises (pf-j8J and (11)-(12). Note that the 



structure of (11 )-(12) for S corresponds purely to (pb; it will not be the covariance 



structure in the posterior inference for /3. 
3.3 Results of Implementation 

For the investigation, we focus on a single LHFI-A model involving the covariates 
DD, log(depth), log(SC) and log(depth)xlog(SC). Bayesian parameter estimates 



^See |Hoff |2QQ9| for a discussion on lack of identifiability. 
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Figure 2: Brooks-Gelman-Rubin diagnostics ( [Spiegelhalter et oL . 2011[ ) for the 



MCMC samples of each matrix element of a block diagonal S. "Sig.b[i^j] 
denotes S^^, i.e. the (i, j')th element of S. Convergence is suggested by a red curve 
approaching 1, and green and blue curves approaching the same constant. 



and credible intervals then are compared to the LHFI-A models from Wu (un 
published! with the same set of covariates. Though, instead of fitting (11)-(|12[) 



which would require an implementation outside of OpenBUGS, we implemented 
a block diagonal S (i.e. ^=0) as a limiting case. Inference for Richibucto latent 
health as a whol^ from this limiting case was very similar to that from a diago- 
nal S as assumed by |Wu ( unpublished[ ), suggesting reasonable robustness of the 



health inference to S. Additionally, the concern of weak identifiability associated 
with a non-diagonal S proves to be somewhat irrelevant for these data, as our two 
independently generated MCMC chains mixed readily after a burn-in of around 
20,000 iterations: parameters of the non-diagonal S required this longer burn- 
in (Fig. |2]), although all other model parameters each required a burn-in of only 
l,00Cp(Fig.|3|). 

Finally, we observe that the significance of the covariates was also virtually 
unaffected by assuming a S that is more complex than that for Wu| ( unpublished) . 



4 An Additional Level in the Latent Regression 

The investigation of Section |3] suggested that the extra complexity from a non- 
trivial dependence structure among metric effects does not help to clarify the re- 



"Health inference as a whole'' here refers to the ranking of sites according to the posterior 
distributions of i^^s. 

^^For a given model, inference for model parameters as a whole was always based on the longest 
burn-in required. 



11 












1000 25000 


50000 75000 1 00000 
iteration 











50000 
iteration 







0.0 0.5 


1 000 25000 


50000 75000 1 00000 
iteration 










1000 25000 


50000 75000 100000 
iteration 






1000 25000 


50000 75000 100000 
iteration 



Figure 3: Trace plots of MCMC iterations (thinned by 100) from the poste- 
rior for fitting (|3])-(|9]); each covariate in ^ has been centred (see Footnote [7|). 
Counter-clockwise from top-left: regression coefficient of the centred interaction 
log(depth) X log(SC), intercept a^, random effect /32, fixed effect latent health 
His, arid standard deviation an- Trace plots for all other non-S model parameters 
show similar patterns that suggest convergence after a burn-in of merely 1,000. 

lationships between latent health and covariates given the current Richibuto data. 
In this section, we revert to the naive independence assumption for /3 but intro- 
duce extra model complexity through an additional level in the regression of latent 
health on covariates. 

Specifically, although the strong correlation between salinity and DD reflects 
ecological reasoning for coastal sea waters entering an estuary, it is the only clear 
empirical relationship detected among the available covariates. Therefore, instead 
of considering salinity and DD to be complementary covariates, we now take salin- 
ity as a response of DD, and in turn, latent health as a response of salinity and the 
remaining covariates from Section 3.3 (Fig.|4]). 

Then, model statements of this section include ([l])-([8]), and additionally. 



^sal,i — <^DD^DD,i + , 

Normal(0, a^) 



CI 2 iid 



(13) 
(14) 



where S = a^I in ds]), I is the identity matrix, and Xi in ^ denotes the vector of 
centred covariates for site i including salinity Xsau (and possibly other covariates) 
but excluding DD XDD,i- Hence, ^ and (p3|) can be collapsed as follows: 



Hi = ao + a'_g^ia3_sal,i + 0^sal<^DD^DD,i + C^saA + £i 



(15) 



where a_sai is cx with o^sai removed, and similarly for x^sai^i- Thus, ([15]) regards 
salinity as an implicit covariate, so that when latent health is explicitly regressed 
on X-sai.i and DD, the implicit covariate decomposes the total error variation into 
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Figure 4: Graphical depiction of regressing salinity on DD as an additional level 
in the hierarchical latent health model, while, as usual, health is the response of 
salinity and other covariates, and AMBI metrics are responses of latent health. 



Var(asai^^i + £i\as^h crh ^h) = ^lai^s + ^h- Hence, a smaller ratio cr|^/(o^sai^^ + ^h) 
reflects a higher contribution from the implicit covariate towards explaining the 
total error variation of the latent health regression. 

We employ the same inverse-Gamma prior for aj^ as well as a^. Univariate 
N(0, 100) priors are employed for ag, and vector components of a, with one 
exception: Cor(asab <^dd) = p 7^ is additionally considered, where p ~ Unif(— 1, 1) 
a priori. 



4.1 Results 

Inference summaries appear in Table [3| in which Models (l)-(3) each comprises 
two levels of covariates (expressions Q-(|8]) and (p3])-(p^), and Models (4) and 
(5) — provided as a comparison — each comprises a single level of covariates 
(expressions ([l])-([8]) only). Posterior means for latent health along with their 95% 
posterior credible intervals (CIs) appear Fig. |5| those for ao, /3j(s)/ ^5/ ^if/ arid a/3 
appear in Fig. |6} 

For the assessment and prediction of latent health, the model parameters of 
main interest include Hi, an, and a, while the rest of the parameters are regarded 
as nuisance ( Wu[ unpublished[ ). 



It is evident from the 95% CIs in Table [3] that, when DD is considered a driver 
of salinity, the two are simultaneously relevant to explaining latent health. In fact. 
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Table 3: Selected summary statistics of posterior draws. Boldfaced CI limits sug- 
gest that the corresponding parameter differs from with high credibility 



95% CI 



Model 



Mean Median 2.5% 97.5% DIC 



(1) sal-on-DD only; 
Cor(asai, Q^dd)=0 



(2) sal-on-DD only; 

Cor(asal,Q^DD)=P 



(3) 



log(depth), log(SC), 
log(depth) X log(SC), 
and sal-on-DD; 
Cor(a)=0 



(4) 



DD only 



(5) 



sal only 



crs 
cth 
O2 



4 



q;dd 
P 

O2 



«0 
Q^dep 
ttsal 

asc 

«DD 
^depxSC 

^2 

q;o 
q;dd 

02 
ttsal 
O2 



-1.56 
0.39 
0.77 
1.12 
0.70 
0.67 
2.09 

0.54 

-1.57 
0.49 
0.59 

-0.94 
1.12 
0.74 
0.68 
2.10 

0.59 

-1.37 
0.16 
0.42 

-0.83 
0.77 
1.88 
1.12 
0.70 
0.59 
2.09 

0.63 

-1.57 
0.37 
1.48 
0.63 
2.09 

-1.57 
0.39 
1.48 
0.67 
2.10 



-1.57 
0.39 
0.77 
1.01 
0.68 
0.65 
2.10 



-1.58 
0.49 
0.59 

-0.96 
1.01 
0.72 
0.66 
2.10 

0.59 

-1.37 
0.16 
0.42 

-0.83 
0.77 
1.88 
1.01 
0.68 
0.57 
2.10 

0.63 

-1.58 
0.37 
1.03 
0.61 
2.10 

-1.57 
0.39 
1.02 
0.66 
2.10 



-3.25 
0.13 
0.54 

0.59 
0.51 
0.48 
-0.11 



0.55 0.31 



-3.29 
0.28 
0.35 

-1.00 

0.59 
0.53 
0.49 
-0.10 

0.37 

-3.08 

-0.63 
0.17 

-1.74 
0.54 
0.28 
0.59 
0.51 
0.41 

-0.10 

0.37 

-3.26 
0.16 
0.35 
0.45 

-0.12 

-3.26 
0.13 

0.35 
0.48 
-0.11 



0.16 
0.64 
1.00 

2.31 
0.98 
0.95 
4.24 

0.76 

0.18 
0.71 
0.79 
-0.78 

2.30 
1.07 
0.96 
4.27 

0.79 

0.36 
0.94 
0.67 
0.08 
1.00 
3.49 
2.29 
0.98 
0.87 
4.24 

0.87 

0.16 
0.58 
5.32 
0.89 
4.23 

0.16 
0.64 

5.32 
0.95 
4.24 



4417 



4419 



4417 



4380 



4379 



is excluded from the 99% CI (not shown) for both agai and q;dd in each of Models 
(l)-(3), suggesting very high credibility for the two-level structure. The CIs from 
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Figure 5: LHFI scores (posterior means marked by and 95% CIs of site health 
(arrows from dark to light), based on Models (4), (5), (1), (2), and (3) of Table [3| 
respectively. |Lu et al, ( 2008| ) partition Richibucto sites into six groups according 
to their benthic community composition: red (lower channel), violet (estuarine 
mouth), green, blue (lower shallow), turqoise (upper shallow), and pink (upper 
channel). 



Model (3) suggest that the interaction log(depth) xlog(SC) is an additional credible 
driver of latent health, complementing the explanatory capacity of salinity-on-DD. 
This is the first time that a scientifically sensible model has been successfully con- 
structed to rigorously identify the collective explanatory capacity of salinity, DD, 
depth, and SC — all regarded a priori as qualitatively important — towards site 
health in the Richibucto ecosystem. 

For Models (l)-(3), the posterior mean for the ratio cr|f/ (<^sai^? + ^1/) I'^rig^s 
from around 0.55 to 0.65, respectively; corresponding 95% CIs suggest a smallest 
ratio for Model (1). Thus, despite the high credibility of the correlation between a^^i 
and iri Model (2) and of the influence on health from (the interaction between) 
depth and SC in Model (3), the least complex Model (1) provides slightly clearer 
evidence for the explanatory capacity of the salinity-on-DD structure. In terms of 
the model's predictive power, the least and most complex among the three mod- 
els share the same deviance information criterion (DIC) ( Spiegelhalter et al\\2.QQ2) 
which is slightly smaller (better) than that of Model (2). This predictive power 
corresponds to observed AMBI metrics (not latent health) being predicted by the 
model. To assess the model's predictive ability for latent health, one could conduct 
a simulation study in which unobservable Hi values are generated then estimated. 



although such an approach for hierarchical models has its shortcomings (Marshall 
& Spiegelhalter[ 2007| ) or requires intensive computations (Dey et ah , 199S) that 



may be impractical. 

Instead, we compare CIs for Hi among models; in Fig. |5| they appear nearly 
identical across all Models (l)-(5), i.e. the inference for health is essentially equally 
credible across various models. Within model, the ranking of sites according to 
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(f) 



Figure 6: Posterior means and 95% CIs for (a) a^, (b) 02, (c) elements of /3^, (d) 
elements of /3_, (e) an, and (f) (plotted on log-scale). 



their LHFI scores and associated CIs do not coincide with the clustering by |Lu| 
et ah (2008 ) based on similarity in benthic community composition (highly corre- 
lated with site location). This indicates that the LHFI approach does not merely 
represent community composition or site location; instead, it rigorously and com- 
prehensively models biotic indicators, abiotic drivers, the abstract notion of health, 
and the relationship among them. Note that the health CIs from the 18 sites mu- 
tually overlap, suggesting that the small dataset does not allow us to distinguish 
sites at a high credible level based on health; this was also the case for those mod- 
els by Wu| ( [unpublished) , all with single-level covariates. Despite (i) suboptimal 
distinguishability and (ii) weaker predictive power for AMBI metrics compared to 
the single-covariate Models (4)-(5), our two-level-covariate Models (l)-(3) clearly 
resolved the earlier counterintuitive phenomenon of covariates not being simul- 
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taneously significant. Indeed, (i) is an improvement over conventional methods 
in quantitative rigour due to the integrated manner from which our uncertainty 
estimates are obtained. Moreover, (ii) is of secondary concern when the response 
of key interest is Hi instead of the metrics. Technical note: The LHFI framework 
is built on the fundamental principles of analysis-of-covariance, so that one can 
only interpret Hi values in a relative sense. However, Chiu et ah ( 2011| ) explain that 
including in the study any site that is qualitatively pre-identified as very healthy 
or very unhealthy would facilitate the interpretation of the magnitude of Hi for 
an individual i. This is slightly different from the approach of [Lopez & Fennessy 



( |2002[ ) who include sites that span the spectrum of individual covariates. 

Finally, aside from nearly identifical CIs for Hi across models. Fig. |6] indicates 
that the five models perform equally well with respect to the uncertainty (width of 
CIs) of various nuisance parameters, but with one exception: two-level-covariate 
models have clearly less uncertainty in their inference for Gj^ (Fig. 6(f) | ). As this 
parameter directly contributes to the uncertainty in the linear predictor v of AMBI 
metrics, Fig. |6(f)| indicates that having two levels of covariates lead to more reliable 
inference for the model as a whole. 



5 Conclusions 

Unlike conventional multimetric health indices, the integrated LHFI approach yields 
health scores, asseses the influence of health drivers, and provides their associated 
uncertainty, all in a single, unified analysis for a given model. 

LHFI models by |Wu| ( |unpublished| ) with single-level covariates were satisfac- 



tory as preliminary models for the 18 Richibucto sites, but lacking the important 
ability to rigorously identify relationships between health and abiotic drivers that 



are deemed ecologically important for the Richibucto estuarine system. Wu (un 



published! proposed, without implementation, two ways to address this issue: (a) 



to introduce a covariance structure on the random metric effects, and (b) to intro- 
duce additional layers of regression given preconceived relationships among the 
covariates. In this paper, we implemented (a) and (b) with AMBI biotic metrics 
only, but the approach would be applicable in principle to combining AMBI and 
ITI biotic metrics. Though, with merely 18 sites in Richibucto, combined AMBI- 
ITI models by |Wu ( |unpublished[ ) suggest that ITI metrics potentially weaken any 



signal in the health-covariate relationship. 

In this paper, we have found (a) to be ineffective. Various block diagonal co- 
variance structures were considered, the most general of which was reported in 
this paper. However, none substantially affected the health inference nor improved 
credible levels of the covariates. 

On the other hand, our efforts on (b) proved to be well rewarded. An addi- 
tional layer of covariates based upon the empirical relationship between salinity 
and distance downstream allowed the model to identify the simultaneous signifi- 
cance of distance and those abiotic covariates that were previously shown by |Wu 



( |unpublished| ) to be significant only when distance was excluded. Moreover, we 
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also show that model inference is more reliable overall when compared to single- 
level-covariate models. Thus, our two-level-covariate modelling framework more 
comprehensively exploits the ecological relationship among health, biotic metrics, 
and abiotic covariates, and it yields less uncertainty in model inference. We im- 
plemented three variants of the two-level-covariate model: (a) salinity-on-distance 
alone, with a priori independent regression coefficients and metric effects, (b) same 
as (a) but assuming bivariate regression coefficients, and (c) same as (a) but includ- 
ing channel depth and silt-clay content (both on the log scale), as well as their in- 
teraction. Overall, (a)-(c) are almost equal in statistical performance, with slightly 
better predictive power of biotic metrics by (a) and (c). Finally, (a) corresponds 
to marginally stronger evidence for the two-level structure between salinity and 
distance to influence site health. 

Although field data, especially biotic data, are costly to collect and process for 
the use in quantitative assessment of ecosystem health, our work has shedded light 
on one practical concern: more than 18 sites and/ or more precise measurements 
on abiotic covariates are needed in order for the LHFI framework to rigorously 
distinguish Richibucto sites according to AMBI and /or ITI metrics as indicators of 
ecosystem health. This point will be considered as part of future health monitoring 
and conservation efforts for the Richibucto estuarine system. 
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